// STL
#include <fstream>                            // std::ofstream
#include <memory>                             // std::unique_ptr<>
#include <utility>                            // std::move()
#include <string>                             // std::string

// Boost
#include <boost/archive/text_oarchive.hpp>    // boost::archive::text_oarchive
#include <boost/serialization/unique_ptr.hpp> // serialize std::unique_ptr<>

// jScience
#include "jScience/linalg.hpp"             // Matrix<>, Vector<>

// stats++
#include "statsxx/machine_learning/Gaussian_process.hpp" // Gaussian process stuff


//
// DESC: "Trains" a GP.
//
std::unique_ptr<Gaussian_process::GaussianProcess> train_GP(
                                                            std::unique_ptr<Gaussian_process::GaussianProcess>  gp,
                                                            // -----
                                                            double                                              sigma_n,
                                                            // -----
                                                            const bool                                          optimize_hyperparameters,
                                                            //------
                                                            const Matrix<double>                               &X_in,
                                                            const Vector<double>                               &x_out,
                                                            // -----
                                                            const std::string                                   gp_file
                                                            )
{
    // NOTE: Return std::unique_ptr<Gaussian_process::GaussianProcess> gp.

    //=========================================================

    //=========================================================
    // INITIALIZATION
    //=========================================================

    std::cout << "here 1" << std::endl;

    if( optimize_hyperparameters )
    {
        std::cout << "here 3x" << std::endl;

        std::vector<double> theta(4);
        // NOTE: sigma_n passed in.
        std::tie(
                 theta,
                 sigma_n
                 ) = Gaussian_process::optimize_hyperparameters(
                                                                gp->k,
                                                                gp->k->theta,
                                                                X_in,
                                                                x_out,
                                                                X_in,
                                                                x_out,
                                                                sigma_n
                                                                );

        std::cout << "here 2x" << std::endl;

        gp->k->assign_hyperparameters(
                                      theta
                                      );

        gp = std::unique_ptr<Gaussian_process::GaussianProcess>(new Gaussian_process::GaussianProcess(gp->k));
    }

    //=========================================================
    // "TRAIN"
    //=========================================================

    gp->add_data(
                 X_in,
                 x_out,
                 sigma_n
                 );

    //=========================================================
    // OUTPUT
    //=========================================================
/*
    // SAVE GP TO ARCHIVE
    {
        std::ofstream ofs(gp_file, std::ios::out);

        boost::archive::text_oarchive oa(ofs);

        oa << gp;
    }
*/
    return std::move(gp);
}
